# spec_analysis.R
# Specification analysis of the Levant VAR models
#
# Patrick Brandt
# pbrandt@utdallas.edu
# Original version: August 2004
# Additional comments: July 2005

# This file evaluates the posterior for the BVAR model for difference
# values of the BVAR hyperparameters.  It reads these hyperparameters
# values from a file, priors.csv and then fits a BVAR and computes the
# posterior measures for each model.


# Set up functions and libraries
library(foreign)
library(MSBVAR)

data <- read.dta("levant.weekly.79-03.dta") 
attach(data)

# Set up KEDS data
KEDS.data <- ts(cbind(a2i,a2p,i2a,p2a,i2p,p2i),
                start=c(1979,15),
                freq=52,
                names=c("A2I","A2P","I2A","P2A","I2P","P2I"))

KEDS <- window(KEDS.data, end=c(1988,50))

################################################################################
# Specification analysis for Lag Length
################################################################################
print(var.lag.specification(KEDS))
print(var.lag.specification(KEDS[,5:6]))

################################################################################
# Specification analysis for the prior
################################################################################
# Build a matrix to store the lag lengths and the prior values with
# the log posteriors.

# Set constants and read in the matrix of priors.
p.max <- 20
p.min <- 4
priors <- read.csv("priors.csv")
k <- 0

results <- matrix(0,nrow=((p.max - p.min + 1)*nrow(priors)), ncol=12)
colnames(results) <- c(colnames(priors[1,1:7]),"Lag","ML","POST","BIC","AIC")


for (p in p.min:p.max)
   { day <- 31-p
     y <- window(KEDS, start=c(1979,day), end=c(1988,50))

     for(i in 1:nrow(priors))
       { prior.tmp <- data.matrix(priors[i,1:7])
         prior.dist <- priors$dist[i]
 
     # Estimate the BVAR for the respective prior
         tmp <- szbvar(y, p, z=NULL, lambda0 = prior.tmp[1],
                       lambda1 = prior.tmp[2],
                       lambda3 = prior.tmp[3],
                       lambda4 = prior.tmp[4],
                       lambda5 = 0,
                       mu5 = 0,
                       mu6 = 0, nu = ncol(y),
                       qm=4,
                       prior=prior.dist,
                       posterior.fit=T)
         # Compute the posterior summaries
         marg.llf <- tmp$marg.llf        
         bic <- 2*marg.llf + prod(dim(tmp$Bhat))*(nrow(y)-p)
         aic <- -2*marg.llf + 2*prod(dim(tmp$Bhat))
         marg.post <- tmp$marg.post
         k <- k+1
         cat("KEDS Data ", k, "\n")
         results[k,] <- c(prior.tmp, p, marg.llf, marg.post, bic,
                          aic)

       }
   }


# Now look at the ordered results, based on the marginal posteriors

best.results <- results[order(results[,10], decreasing=T),]
print(best.results)

# Clean up environment and save results
detach(data)
save.image(file="Levant-spec-analysis.RData", compress=T)


# End of file
